Setup

knitr::opts_chunk$set(echo = TRUE)

library(tidyverse)
library(R.utils)
library(wbCorr)
library(readxl)
library(kableExtra)
library(brms)
library(bayesplot)
library(beepr)


source(file.path('Functions', 'ReportModels.R'))
source(file.path('Functions', 'PrettyTables.R'))
source(file.path('Functions', 'ReportMeasures.R'))
source(file.path('Functions', 'PrepareData.R'))
system("shutdown /a")
## [1] 1116
# Set options for analysis
use_mi = FALSE
shutdown = FALSE
report_ordinal = FALSE

options(
  dplyr.print_max = 100, 
  brms.backend = 'cmdstan',
  brms.file_refit = ifelse(use_mi, 'never', 'on_change'),
  #brms.file_refit = 'always',
  error = function() beepr::beep(sound = 5)
)
df <- openxlsx::read.xlsx('long.xlsx')
df_original <- df

df_list <- prepare_data(df, use_mi = use_mi)

Constructing scales Re-coding pusing reshaping data (4field) centering data within and between

df_double <- df_list[[1]]
df_full <- df_list[[2]]

Descriptives

Sample Statistics

# Income (mean of both partner's report)
merge_income <- function(income1, income2) {
  merged_income <- numeric(length(income1))
  
  # Loop through each pair of incomes
  for (i in seq_along(income1)) {
    # Handle NA
    if (is.na(income1[i])) {
      merged_income[i] <- income2[i]
    } 
    
    else if (is.na(income2[i])) {
      merged_income[i] <- income1[i]
    }
    
    # if both are informative, take mean and round
    else if (income1[i] %in% 1:6 && income2[i] %in% 1:6) {
      merged_income[i] <- round((income1[i] + income2[i]) / 2)
    }
    
    # if one is informative and the other not, use the informative one. 
    else if (income1[i] %in% 1:6) {
      merged_income[i] <- income1[i]
    }
    else if (income2[i] %in% 1:6) {
      merged_income[i] <- income2[i]
    }
    
    # Now we only have cases, where both are either 70 or 99. We simply report "undisclosed". 
    else {
      merged_income[i] <- 99
    }
  }
  
  # Convert to factor
  merged_income <- factor(
    merged_income, levels = c(1,2,3,4,5,6,70,99), 
    labels = c(
      "up to CHF 2'000.-", 
      "CHF 2'001.- to CHF 4'000.-",
      "CHF 4'001.- to CHF 6'000.-",
      "CHF 6'001.- to CHF 8'000.-",
      "CHF 8'001.- to CHF 10'000.-", 
      "above CHF 10'000.-", 
      "I don't know",
      "Undisclosed"
  )
  )
  
  return(merged_income)
}



df_sample_report <- df_full %>%
  group_by(coupleID) %>%
  arrange(userID) %>%
  # Computing couple level variables
  mutate(
    Household_Income = merge_income(first(pre_income_1), last(pre_income_1)),
    
    reldur = pre_rel_duration_m / 12 + pre_rel_duration_y,
    Relationship_duration = mean(reldur, na.rm = TRUE),
    
    habdur = pre_hab_duration_m / 12 + pre_hab_duration_y,
    Cohabiting_duration = mean(habdur, na.rm = TRUE),
    
    Marital_status = factor(
      case_when(
        all(pre_mat_stat == 1) ~ "Married",
        any(pre_mat_stat == 1) ~ "One Partner Married",
        TRUE ~ "Not Married"
      )
    ),
    
    Have_children = factor(
      (first(pre_child_option) + last(pre_child_option)) > 0, 
      levels = c(FALSE, TRUE), 
      labels = c('Have Children', 'No Children')),
    
    Gender = factor(
      gender, 
      levels = c(1,2,3), 
      labels = c('Male','Female', 'Other')),
    

    Couple_type = as.factor(
      case_when(
        first(Gender) == last(Gender) & first(Gender) == 'Male' ~ 'Same-Sex Couple (Male)',
        first(Gender) == last(Gender) & first(Gender) == 'Female' ~ 'Same-Sex Couple (Female)',
        TRUE ~ 'Mixed-sex Couple'
      )
    )
  ) %>%
  ungroup() %>%
  # Individual level variables
  mutate(
    Age = pre_age,
    Handedness = factor(
      pre_handedness, 
      levels = c(0, 1, 2), 
      c('Right','Left', 'Ambidextrous')),
  Highest_Education = factor(
    pre_education, 
    levels = c(1,2,3,4,5,6,7), 
    labels = c(
      "(still) no school diploma",
      "compulsory education (9 years)",
      "vocational training (apprenticeship)",
      "Matura (university entrance qualification)",
      "Bachelor's degree", 
      "Master's degree",
      "Doctorate degree"
      )
    ),
  BMI = pre_weight / ((pre_height / 100)^2) # to meters
  ) %>%
  select(c(Relationship_duration, Cohabiting_duration, Couple_type, Household_Income, 
            Marital_status, Have_children, 
            Gender, Age, Handedness, Highest_Education, BMI))


sample_table <- report_measures(df_sample_report, ICC = F)
sample_table$n_Obs <- as.numeric(sample_table$n_Obs) / 55
rownames(sample_table) <- NULL

n_couple_vars <- 17
sample_table$n_Obs[1:n_couple_vars] <- sample_table$n_Obs[1:n_couple_vars] / 2

packing_sample <- list(
    "Couple level variables (38 couples)" = 
      c(1, n_couple_vars),
    "Individual level variables (76 individuals)" 
    = c(n_couple_vars+1, nrow(sample_table))
    ) 

df_sample_summary <- print_df(
  sample_table,
  rows_to_pack = packing_sample
)

export_xlsx(df_sample_summary,
            file.path('Output', 'SampleDescription.xlsx'),
            merge_option = 'both',
            rows_to_pack = packing_sample,
            colwidths = c(20,35,7,7,7,7,7,10)
            )

df_sample_summary
Variable Level n_Obs percentage_Obs Missing Mean SD Range
Couple level variables (38 couples)
Relationship_duration NA 38 NA 0% 9.23 9.03 0.58-36.00
Cohabiting_duration NA 38 NA 0% 7.53 9.14 0.25-33.00
Couple_type Same-Sex Couple (Male) 1 3% NA NA NA NA
Couple_type Same-Sex Couple (Female) 1 3% NA NA NA NA
Couple_type Mixed-sex Couple 36 95% NA NA NA NA
Household_Income I don’t know 0 0% NA NA NA NA
Household_Income up to CHF 2’000.- 2 5% NA NA NA NA
Household_Income CHF 2’001.- to CHF 4’000.- 3 8% NA NA NA NA
Household_Income CHF 6’001.- to CHF 8’000.- 3 8% NA NA NA NA
Household_Income CHF 4’001.- to CHF 6’000.- 5 13% NA NA NA NA
Household_Income Undisclosed 5 13% NA NA NA NA
Household_Income CHF 8’001.- to CHF 10’000.- 8 21% NA NA NA NA
Household_Income above CHF 10’000.- 12 32% NA NA NA NA
Marital_status Married 13 34% NA NA NA NA
Marital_status Not Married 25 66% NA NA NA NA
Have_children No Children 11 29% NA NA NA NA
Have_children Have Children 27 71% NA NA NA NA
Individual level variables (76 individuals)
Gender Other 0 0% NA NA NA NA
Gender Male 38 50% NA NA NA NA
Gender Female 38 50% NA NA NA NA
Age NA 76 NA 0% 34.01 10.96 19.00-60.00
Handedness Ambidextrous 1 1% NA NA NA NA
Handedness Left 11 14% NA NA NA NA
Handedness Right 64 84% NA NA NA NA
Highest_Education (still) no school diploma 0 0% NA NA NA NA
Highest_Education compulsory education (9 years) 0 0% NA NA NA NA
Highest_Education Doctorate degree 0 0% NA NA NA NA
Highest_Education vocational training (apprenticeship) 8 11% NA NA NA NA
Highest_Education Master’s degree 21 28% NA NA NA NA
Highest_Education Bachelor’s degree 23 30% NA NA NA NA
Highest_Education Matura (university entrance qualification) 24 32% NA NA NA NA
BMI NA 76 NA 0% 24.94 4.08 16.37-33.95

Measures

Focal Variables

main_constructs <- c("persuasion", "pressure","pushing", 
                     "pa_sub", "pa_obj", "aff", "reactance"
                     )

main_descriptives <- report_measures(
  data = df_full, 
  measures = main_constructs,
  ICC = TRUE, 
  cluster_var = df_full$userID)

openxlsx::write.xlsx(
  main_descriptives, 
  file.path('Output', 'DescriptivesMain.xlsx')
  )

print_df(main_descriptives)
Variable n_Obs Missing Mean SD Range ICC
persuasion 4180 6% 0.42 1.08 0.00- 5.00 0.23
pressure 4180 6% 0.12 0.62 0.00- 5.00 0.58
pushing 4180 6% 0.16 0.66 0.00- 5.00 0.11
pa_sub 4180 6% 30.40 54.78 0.00-720.00 0.16
pa_obj 4180 11% 144.41 117.81 5.75-971.25 0.18
aff 4180 6% 4.83 1.14 1.00- 6.00 0.41
reactance 4180 81% 0.79 1.32 0.00- 5.00 0.47

All Variables

all_constructs <- c(
  main_constructs,
  "day",
  "weartime",
  "isWeekend",
  "plan",
  "studyGroup",
  "support",
  "got_JITAI",
  "skilled_support"
)


all_descriptives <- report_measures(df_full, all_constructs, ICC = F)

openxlsx::write.xlsx(
  all_descriptives, 
  file.path('Output', 'DescriptivesAll.xlsx')
)

print_df(all_descriptives)
Variable Level n_Obs percentage_Obs Missing Mean SD Range
persuasion NA 4180 NA 6% 0.42 1.08 0.00- 5.00
pressure NA 4180 NA 6% 0.12 0.62 0.00- 5.00
pushing NA 4180 NA 6% 0.16 0.66 0.00- 5.00
pa_sub NA 4180 NA 6% 30.40 54.78 0.00- 720.00
pa_obj NA 4180 NA 11% 144.41 117.81 5.75- 971.25
aff NA 4180 NA 6% 4.83 1.14 1.00- 6.00
reactance NA 4180 NA 81% 0.79 1.32 0.00- 5.00
day NA 4180 NA 0% 0.50 0.29 0.00- 1.00
weartime NA 4180 NA 1% 1057.42 384.29 0.00-1440.00
isWeekend Weekend 1216 29% NA NA NA NA
isWeekend Weekday 2964 71% NA NA NA NA
plan missing 238 6% NA NA NA NA
plan Plan 1860 44% NA NA NA NA
plan No plan 2082 50% NA NA NA NA
studyGroup last 3 weeks interventions 1320 32% NA NA NA NA
studyGroup Allways inerventions 1430 34% NA NA NA NA
studyGroup First 3 weeks interventions 1430 34% NA NA NA NA
support NA 4180 NA 6% 0.91 1.52 0.00- 5.00
got_JITAI JITAI received 585 14% NA NA NA NA
got_JITAI No JITAI 3595 86% NA NA NA NA
skilled_support Days before Intervention 1036 25% NA NA NA NA
skilled_support Days after Intervention 3144 75% NA NA NA NA

Correlations

cors <- wbCorr(df_full[,c(main_constructs)], df_full$coupleID, method = 'spearman')

main_cors <- summary(cors, 'wb')$merged_wb


openxlsx::write.xlsx(
  main_cors, 
  file.path('Output', 'Correlations.xlsx')
)

print_df(main_cors, width = '7em')
persuasion pressure pushing pa_sub pa_obj aff reactance
persuasion [0.20] 0.21*** 0.40*** 0.17*** 0.07*** 0.00 -0.05
pressure 0.30 [0.55] 0.28*** -0.03 -0.04* -0.01 0.26***
pushing 0.63*** 0.40* [0.08] 0.14*** 0.05** 0.01 0.07*
pa_sub 0.27 -0.10 0.24 [0.15] 0.31*** 0.20*** -0.04
pa_obj 0.13 -0.08 0.27 0.51*** [0.13] 0.19*** 0.06
aff 0.28 -0.07 0.29 0.52*** 0.22 [0.28] -0.01
reactance 0.18 0.23 0.11 0.06 0.38* -0.12 [0.42]

Within-person correlations are above the diagonal and between-person correlations are below the diagonal.
On the diagonal are intraclass correlations (ICCs)

Modelling

# For indistinguishable Dyads
model_rows_fixed <- c(
    'Intercept', 
    # '-- WITHIN PERSON MAIN EFFECTS --', 
    'persuasion_self_cw', 
    'persuasion_partner_cw', 
    'pressure_self_cw', 
    'pressure_partner_cw', 
    'pushing_self_cw', 
    'pushing_partner_cw', 
    'day', 
    'weartime_self_cw',
    
    # '-- BETWEEN PERSON MAIN EFFECTS',
    'persuasion_self_cb',
    'persuasion_partner_cb',
    'pressure_self_cb',
    'pressure_partner_cb',
    'pushing_self_cb',
    'pushing_partner_cb',
    'weartime_self_cb'
  )


model_rows_fixed_ordinal <- c(
  model_rows_fixed[1],
  'Intercept[1]',
  'Intercept[2]',
  'Intercept[3]',
  'Intercept[4]',
  'Intercept[5]',
  model_rows_fixed[2:length(model_rows_fixed)]
)

model_rows_random <- c(
  # '--------------',
  # '-- RANDOM EFFECTS --',
  'sd(Intercept)', 
  'sd(persuasion_self_cw)',
  'sd(persuasion_partner_cw)',
  'sd(pressure_self_cw)',
  'sd(pressure_partner_cw)',
  'sd(pushing_self_cw)',
  'sd(pushing_partner_cw)',
  # '-- CORRELATION STRUCTURE -- ', 
  'ar[1]', 
  'nu',
  'shape',
  'sderr',
  'sigma'
)

model_rows_random_ordinal <- c(model_rows_random,'disc')
# For indistinguishable Dyads
model_rownames_fixed <- c(
    'Intercept', 
    # '-- WITHIN PERSON MAIN EFFECTS --', 
    'Daily received persuasion target -> target', 
    'Daily received persuasion target -> agent', 
    'Daily received pressure target -> target', 
    'Daily received pressure target -> agent', 
    'Daily received pushing target -> target', 
    'Daily received pushing target -> agent', 
    'Day', 
    'Daily weartime',
    
    # '-- BETWEEN PERSON MAIN EFFECTS',
    'Mean received persuasion target -> target',
    'Mean received persuasion target -> agent',
    'Mean received pressure target -> target',
    'Mean received pressure target -> agent',
    'Mean received pushing target -> target',
    'Mean received pushing target -> agent',
    'Mean weartime'
  )


model_rownames_fixed_ordinal <- c(
  model_rownames_fixed[1],
  'Intercept[1]',
  'Intercept[2]',
  'Intercept[3]',
  'Intercept[4]',
  'Intercept[5]',
  model_rownames_fixed[2:length(model_rownames_fixed)]
)

model_rownames_random <- c(
  # '--------------',
  # '-- RANDOM EFFECTS --',
  'sd(Intercept)', 
  'sd(Daily received persuasion target -> target)', 
  'sd(Daily received persuasion target -> agent)', 
  'sd(Daily received pressure target -> target)', 
  'sd(Daily received pressure target -> agent)', 
  'sd(Daily received pushing target -> target)', 
  'sd(Daily received pushing target -> agent)', 
  # '-- CORRELATION STRUCTURE -- ', 
  'ar[1]', 
  'nu',
  'shape',
  'sderr',
  'sigma'
)

model_rownames_random_ordinal <- c(model_rownames_random,'disc')
# For indistinguishable Dyads
model_rownames_fixed <- c(
    "Intercept", 
    # "-- WITHIN PERSON MAIN EFFECTS --", 
    "Daily persuasion experienced", 
    "Daily persuasion utilized (partner's view)", # OR partner received
    "Daily pressure experienced", 
    "Daily pressure utilized (partner's view)", 
    "Daily pushing experienced", 
    "Daily pushing utilized (partner's view)", 
    "Day", 
    "Daily weartime",
    
    # "-- BETWEEN PERSON MAIN EFFECTS",
    "Mean persuasion experienced", 
    "Mean persuasion utilized (partner's view)", 
    "Mean pressure experienced", 
    "Mean pressure utilized (partner's view)", 
    "Mean pushing experienced", 
    "Mean pushing utilized (partner's view)", 
    "Mean weartime"
  )


model_rownames_fixed_ordinal <- c(
  model_rownames_fixed[1],
  'Intercept[1]',
  'Intercept[2]',
  'Intercept[3]',
  'Intercept[4]',
  'Intercept[5]',
  model_rownames_fixed[2:length(model_rownames_fixed)]
)

model_rownames_random <- c(
  # '--------------',
  # '-- RANDOM EFFECTS --',
  'sd(Intercept)', 
  "sd(Daily persuasion experienced)", 
  "sd(Daily persuasion utilized (partner's view))", # OR partner received
  "sd(Daily pressure experienced)", 
  "sd(Daily pressure utilized (partner's view))", 
  "sd(Daily pushing experienced)", 
  "sd(Daily pushing utilized (partner's view))", 
  # '-- CORRELATION STRUCTURE -- ', 
  'ar[1]', 
  'nu',
  'shape',
  'sderr',
  'sigma'
)

model_rownames_random_ordinal <- c(model_rownames_random,'disc')
rows_to_pack <- list(
  "Within-Person Effects" = c(2,9),
  "Between-Person Effects" = c(10,16),
  "Random Effects" = c(17, 23), 
  "Additional Parameters" = c(24,28)
  )


rows_to_pack_ordinal <- list(
  "Intercepts" = c(1,6),
  "Within-Person Effects" = c(2+5,9+5),
  "Between-Person Effects" = c(10+5,16+5),
  "Random Effects" = c(17+5, 23+5), 
  "Additional Parameters" = c(24+5,28+6)
  )

Subjective MVPA

range(df_double$pa_sub, na.rm = T) 
## [1]   0 720
hist(df_double$pa_sub, breaks = 100)

Modelling using the gaussian family fails. Due to the many zeros, transformations won’t help estimating the models. We employ the negative binomial family.

formula <- bf(
  pa_sub ~ 
    persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw +
    
    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +
    pushing_self_cb + pushing_partner_cb +
    
    day + 
    
    # Random effects
    (persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw | coupleID),

  autocor = ~ ar(time = day, gr = coupleID:userID, p = 1)
)




prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b"),
  brms::set_prior("normal(0, 50)", class = "Intercept", lb = 0),
  
  brms::set_prior("normal(0, 10)", class = "sd", group = "coupleID", lb = 0),
  
  brms::set_prior("cauchy(0, 5)", class = "ar", lb = -1, ub = 1),
  brms::set_prior("cauchy(0, 20)", class = "shape"), 
  brms::set_prior("cauchy(0, 10)", class='sderr')
)


#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pa_sub <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = brms::negbinomial(),
  #control = list(adapt_delta = 0.99),
  iter = 5000,
  warmup = 2000,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", "pa_sub")
)
pp_check(pa_sub, type='hist')
## Using 10 posterior draws for ppc type 'hist' by default.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

pp_check(pa_sub)
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

loo(pa_sub)
## 
## Computed from 12000 by 3736 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo -12067.1 177.4
## p_loo        38.8   3.2
## looic     24134.2 354.8
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.4, 2.0]).
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. ESS
## (-Inf, 0.7]   (good)     3731  99.9%   464     
##    (0.7, 1]   (bad)         5   0.1%   <NA>    
##    (1, Inf)   (very bad)    0   0.0%   <NA>    
## See help('pareto-k-diagnostic') for details.
plot(pa_sub, ask = FALSE)

summarize_brms(
  pa_sub, 
  model_rows_fixed = model_rows_fixed,
  model_rows_random = model_rows_random,
  model_rownames_fixed = model_rownames_fixed,
  model_rownames_random = model_rownames_random,
  exponentiate = T) %>%
  print_df(rows_to_pack = rows_to_pack)
IRR l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 26.12* 18.60 36.84 1.003 1947.31 4164.19
Within-Person Effects
Daily persuasion experienced 0.99 0.89 1.10 1.001 15403.95 9933.47
Daily persuasion utilized (partner’s view) 0.93 0.84 1.03 1.000 18116.11 9513.42
Daily pressure experienced 1.02 0.78 1.39 1.000 17207.47 8708.50
Daily pressure utilized (partner’s view) 0.64* 0.51 0.82 1.000 15112.84 8859.28
Daily pushing experienced 1.03 0.89 1.22 1.001 15519.65 9498.99
Daily pushing utilized (partner’s view) 1.12 0.96 1.31 1.000 14250.05 9140.63
Day 1.07 0.78 1.47 1.001 18058.09 8552.03
Daily weartime NA NA NA NA NA NA
Between-Person Effects
Mean persuasion experienced 1.22 0.73 2.08 1.000 6488.18 8273.13
Mean persuasion utilized (partner’s view) 1.25 0.73 2.17 1.001 6272.04 8209.06
Mean pressure experienced 1.79 0.89 3.65 1.000 8785.18 9147.19
Mean pressure utilized (partner’s view) 0.52 0.25 1.12 1.000 7770.75 8493.28
Mean pushing experienced 0.40* 0.16 1.00 1.000 7373.08 8987.58
Mean pushing utilized (partner’s view) 0.57 0.26 1.29 1.001 8170.66 8685.38
Mean weartime NA NA NA NA NA NA
Random Effects
sd(Intercept) 0.67 0.49 0.89 1.00 3614.10 6268.21
sd(Daily persuasion experienced) 0.20 0.05 0.35 1.00 2887.62 2333.56
sd(Daily persuasion utilized (partner’s view)) 0.19 0.04 0.34 1.00 3209.49 2177.31
sd(Daily pressure experienced) 0.16 0.01 0.48 1.00 7612.01 6293.85
sd(Daily pressure utilized (partner’s view)) 0.14 0.01 0.43 1.00 8778.82 6999.62
sd(Daily pushing experienced) 0.24 0.02 0.52 1.00 3387.34 3385.96
sd(Daily pushing utilized (partner’s view)) 0.15 0.01 0.35 1.00 5772.56 5034.70
Additional Parameters
ar[1] 0.02 -0.94 0.94 1.00 14024.56 7628.24
nu NA NA NA NA NA NA
shape 0.14 0.13 0.14 1.00 17219.02 8584.43
sderr 0.05 0.00 0.13 1.00 7232.94 5199.90
sigma NA NA NA NA NA NA

Device Based MVPA

range(df_double$pa_obj, na.rm = T) 
## [1]   5.75 971.25
hist(df_double$pa_obj, breaks = 50)

df_double$pa_obj_log <- log(df_double$pa_obj)

hist(df_double$pa_obj_log, breaks = 50)

We tried negative binomial here as well for consistency, but the model did not converge. Poisson also did not work. As we have no zeros in this distribution, we log transform.

formula <- bf(
  pa_obj_log ~ 
    persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw +
    
    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +
    pushing_self_cb + pushing_partner_cb +
    
    day + weartime_self_cw + weartime_self_cb +
    
    # Random effects
    (persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw | coupleID),

  autocor = ~ ar(time = day, gr = coupleID:userID, p = 1)
)



prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b"),
  brms::set_prior("normal(0, 50)", class = "Intercept", lb = 0),
  
  brms::set_prior("normal(0, 10)", class = "sd", group = "coupleID", lb = 0),
  
  brms::set_prior("cauchy(0, 5)", class = "ar", lb = -1, ub = 1),
  brms::set_prior("cauchy(0, 10)", class = "sigma", lb = 0)
)


#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pa_obj_log <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = gaussian(),
  #control = list(adapt_delta = 0.99),
  iter = 5000,
  warmup = 2000,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", "pa_obj_log")
)
# plotting with the first imputed dataset. 
pp_check(pa_obj_log, type='hist')
## Using 10 posterior draws for ppc type 'hist' by default.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

pp_check(pa_obj_log)
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

loo(pa_obj_log)
## 
## Computed from 12000 by 3337 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo  -2811.5  55.6
## p_loo        99.8   4.7
## looic      5623.0 111.2
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.4, 2.4]).
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. ESS
## (-Inf, 0.7]   (good)     3335  99.9%   721     
##    (0.7, 1]   (bad)         2   0.1%   <NA>    
##    (1, Inf)   (very bad)    0   0.0%   <NA>    
## See help('pareto-k-diagnostic') for details.
plot(pa_obj_log, ask = FALSE)

summarize_brms(
  pa_obj_log, 
  model_rows_fixed = model_rows_fixed,
  model_rows_random = model_rows_random,
  model_rownames_fixed = model_rownames_fixed,
  model_rownames_random = model_rownames_random,
  exponentiate = T) %>%
  print_df(rows_to_pack = rows_to_pack)
exp(Est.) l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 117.56* 106.28 129.91 1.001 3110.81 5862.95
Within-Person Effects
Daily persuasion experienced 1.00 0.97 1.02 1.000 24523.62 9306.48
Daily persuasion utilized (partner’s view) 1.02* 1.00 1.05 1.000 27998.42 8695.83
Daily pressure experienced 0.98 0.93 1.04 1.000 24809.95 8629.74
Daily pressure utilized (partner’s view) 1.01 0.96 1.07 1.000 27658.56 9513.25
Daily pushing experienced 1.01 0.98 1.04 1.000 24621.20 9157.59
Daily pushing utilized (partner’s view) 0.98 0.95 1.01 1.000 22973.28 9805.58
Day 1.00 0.92 1.08 1.000 28118.40 8854.86
Daily weartime 1.00 1.00 1.00 1.001 12950.18 7997.18
Between-Person Effects
Mean persuasion experienced 0.88 0.74 1.04 1.000 7540.04 8894.90
Mean persuasion utilized (partner’s view) 1.07 0.89 1.27 1.000 7041.34 8846.56
Mean pressure experienced 1.12 0.90 1.40 1.000 13044.62 9914.54
Mean pressure utilized (partner’s view) 0.77* 0.62 0.95 1.000 9459.61 9061.46
Mean pushing experienced 1.09 0.82 1.45 1.000 9652.42 10168.43
Mean pushing utilized (partner’s view) 1.23 0.97 1.57 1.000 11575.82 9854.78
Mean weartime 1.00* 1.00 1.00 1.000 13775.51 10920.13
Random Effects
sd(Intercept) 0.27 0.20 0.35 1.00 3865.48 6533.30
sd(Daily persuasion experienced) 0.06 0.03 0.09 1.00 8703.02 8246.50
sd(Daily persuasion utilized (partner’s view)) 0.05 0.02 0.08 1.00 7875.45 7043.60
sd(Daily pressure experienced) 0.06 0.00 0.15 1.00 4891.98 5643.63
sd(Daily pressure utilized (partner’s view)) 0.04 0.00 0.10 1.00 8297.51 6639.51
sd(Daily pushing experienced) 0.08 0.01 0.15 1.00 2976.64 3239.61
sd(Daily pushing utilized (partner’s view)) 0.04 0.00 0.10 1.00 4888.28 6292.34
Additional Parameters
ar[1] 0.28 0.24 0.31 1.00 28353.48 9019.94
nu NA NA NA NA NA NA
shape NA NA NA NA NA NA
sderr NA NA NA NA NA NA
sigma 0.55 0.54 0.57 1.00 23241.62 8371.73

Affect

range(df_double$aff, na.rm = T) 
## [1] 1 6
hist(df_double$aff, breaks = 15)

formula <- bf(
  aff ~ 
    persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw +
    
    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +
    pushing_self_cb + pushing_partner_cb +
    
    day +
    
    # Random effects
    (persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw | coupleID),

  autocor = ~ ar(time = day, gr = coupleID:userID, p = 1)
)



prior1 <- c(
  brms::set_prior("normal(0, 5)", class = "b"),
  brms::set_prior("normal(0, 20)", class = "Intercept", lb=0, ub=6), # range of the outcome scale
  
  brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0),

  brms::set_prior("cauchy(0, 5)", class = "ar", lb = -1, ub = 1),
  brms::set_prior("cauchy(0, 10)", class = "sigma", lb = 0)
  
)

#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

mood <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = gaussian(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  iter = 5000,
  warmup = 2000,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", "mood")
)
pp_check(mood, type='hist')
## Using 10 posterior draws for ppc type 'hist' by default.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

pp_check(mood)
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

loo(mood)
## 
## Computed from 12000 by 3736 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo  -4819.4  63.7
## p_loo        91.8   4.3
## looic      9638.9 127.4
## ------
## MCSE of elpd_loo is 0.1.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.4, 2.3]).
## 
## All Pareto k estimates are good (k < 0.7).
## See help('pareto-k-diagnostic') for details.
plot(mood, ask = FALSE)

summarize_brms(
  mood, 
  model_rows_fixed = model_rows_fixed,
  model_rows_random = model_rows_random,
  model_rownames_fixed = model_rownames_fixed,
  model_rownames_random = model_rownames_random,
  exponentiate = F) %>%
  print_df(rows_to_pack = rows_to_pack)
b l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 4.92* 4.71 5.13 1.001 1846.44 3916.97
Within-Person Effects
Daily persuasion experienced 0.02 -0.01 0.05 1.001 19881.67 9474.31
Daily persuasion utilized (partner’s view) -0.02 -0.05 0.01 1.000 19796.07 9102.58
Daily pressure experienced -0.03 -0.11 0.04 1.000 22197.01 8870.65
Daily pressure utilized (partner’s view) 0.04 -0.04 0.11 1.001 23180.76 9713.33
Daily pushing experienced -0.05* -0.10 -0.01 1.001 22110.08 9431.94
Daily pushing utilized (partner’s view) 0.01 -0.04 0.05 1.001 20430.93 9638.76
Day -0.08 -0.22 0.07 1.001 23140.20 9261.49
Daily weartime NA NA NA NA NA NA
Between-Person Effects
Mean persuasion experienced 0.27 -0.04 0.59 1.001 6173.74 7859.18
Mean persuasion utilized (partner’s view) 0.08 -0.23 0.40 1.001 6049.30 8671.07
Mean pressure experienced -0.01 -0.37 0.35 1.000 10825.49 9325.16
Mean pressure utilized (partner’s view) -0.11 -0.51 0.27 1.000 9149.64 8575.72
Mean pushing experienced -0.46 -0.98 0.07 1.000 7635.42 8724.00
Mean pushing utilized (partner’s view) -0.19 -0.65 0.26 1.000 8455.22 9263.59
Mean weartime NA NA NA NA NA NA
Random Effects
sd(Intercept) 0.63 0.49 0.81 1.00 2430.18 5247.64
sd(Daily persuasion experienced) 0.03 0.00 0.07 1.00 5189.44 5732.22
sd(Daily persuasion utilized (partner’s view)) 0.05 0.00 0.11 1.00 3054.58 3364.81
sd(Daily pressure experienced) 0.09 0.00 0.26 1.00 4507.06 6043.70
sd(Daily pressure utilized (partner’s view)) 0.15 0.01 0.33 1.00 3492.99 4501.00
sd(Daily pushing experienced) 0.09 0.01 0.16 1.00 4004.65 2354.41
sd(Daily pushing utilized (partner’s view)) 0.07 0.01 0.16 1.00 3618.65 4004.24
Additional Parameters
ar[1] 0.45 0.42 0.48 1.00 17929.11 9470.12
nu NA NA NA NA NA NA
shape NA NA NA NA NA NA
sderr NA NA NA NA NA NA
sigma 0.87 0.85 0.89 1.00 20518.52 8704.26

Reactance

Gaussian

range(df_double$reactance, na.rm = T) 
## [1] 0 5
hist(df_double$reactance, breaks = 6)

formula <- bf(
  reactance ~ 
    persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw +
    
    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +
    pushing_self_cb + pushing_partner_cb +
    
    day +
    
    # Random effects
    (persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw | coupleID),

  autocor = ~ ar(time = day, gr = coupleID:userID, p = 1)
)



prior1 <- c(
  brms::set_prior("normal(0, 5)", class = "b"),
  brms::set_prior("normal(0, 20)", class = "Intercept", lb=0, ub=5), # range of the outcome scale
  
  brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0),

  brms::set_prior("cauchy(0, 5)", class = "ar", lb = -1, ub = 1),
  brms::set_prior("cauchy(0, 10)", class = "sigma", lb = 0)
)




#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

reactance <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = gaussian(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  iter = 5000,
  warmup = 2000,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", "reactance")
)
pp_check(reactance, type='hist')
## Using 10 posterior draws for ppc type 'hist' by default.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

pp_check(reactance)
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

loo(reactance)
## 
## Computed from 12000 by 756 log-likelihood matrix.
## 
##          Estimate   SE
## elpd_loo  -1076.1 33.9
## p_loo        85.3  8.1
## looic      2152.3 67.8
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.5, 1.9]).
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. ESS
## (-Inf, 0.7]   (good)     750   99.2%   180     
##    (0.7, 1]   (bad)        6    0.8%   <NA>    
##    (1, Inf)   (very bad)   0    0.0%   <NA>    
## See help('pareto-k-diagnostic') for details.
plot(reactance, ask = FALSE)

summarize_brms(
  reactance, 
  model_rows_fixed = model_rows_fixed,
  model_rows_random = model_rows_random,
  model_rownames_fixed = model_rownames_fixed,
  model_rownames_random = model_rownames_random,
  exponentiate = F) %>%
  print_df(rows_to_pack = rows_to_pack)
b l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 0.54* 0.27 0.81 1.001 6164.65 7659.61
Within-Person Effects
Daily persuasion experienced 0.03 -0.02 0.08 1.000 23483.84 9703.39
Daily persuasion utilized (partner’s view) 0.00 -0.06 0.06 1.000 22569.49 10403.33
Daily pressure experienced -0.02 -0.12 0.07 1.000 22173.63 9192.27
Daily pressure utilized (partner’s view) 0.00 -0.13 0.12 1.000 21194.04 9735.53
Daily pushing experienced 0.03 -0.03 0.09 1.000 21514.93 9640.07
Daily pushing utilized (partner’s view) -0.02 -0.10 0.06 1.000 22943.85 9395.86
Day -0.15 -0.42 0.12 1.001 22240.29 9710.18
Daily weartime NA NA NA NA NA NA
Between-Person Effects
Mean persuasion experienced -0.12 -0.49 0.26 1.000 11500.61 9822.64
Mean persuasion utilized (partner’s view) 0.16 -0.24 0.56 1.001 12913.02 9958.06
Mean pressure experienced 0.18 -0.25 0.61 1.000 13897.18 9375.82
Mean pressure utilized (partner’s view) -0.09 -0.57 0.40 1.000 13183.70 9640.78
Mean pushing experienced -0.17 -0.76 0.42 1.000 11261.18 10322.17
Mean pushing utilized (partner’s view) -0.11 -0.71 0.52 1.000 12414.44 10033.75
Mean weartime NA NA NA NA NA NA
Random Effects
sd(Intercept) 0.58 0.44 0.77 1.00 6151.55 7785.35
sd(Daily persuasion experienced) 0.06 0.00 0.15 1.00 3308.48 6441.77
sd(Daily persuasion utilized (partner’s view)) 0.04 0.00 0.12 1.00 7255.40 6743.75
sd(Daily pressure experienced) 0.43 0.27 0.66 1.00 8406.83 9383.52
sd(Daily pressure utilized (partner’s view)) 0.21 0.01 0.58 1.00 3467.39 6248.89
sd(Daily pushing experienced) 0.09 0.00 0.22 1.00 3205.32 6605.58
sd(Daily pushing utilized (partner’s view)) 0.04 0.00 0.14 1.00 7598.29 7477.10
Additional Parameters
ar[1] 0.01 -0.07 0.10 1.00 17003.80 9563.44
nu NA NA NA NA NA NA
shape NA NA NA NA NA NA
sderr NA NA NA NA NA NA
sigma 0.94 0.89 1.00 1.00 12850.87 8795.47 
hypothesis(reactance, "pressure_self_cw > pushing_self_cw")
## Hypothesis Tests for class b:
##                 Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio Post.Prob Star
## 1 (pressure_self_cw... > 0    -0.06      0.06    -0.16     0.04       0.23      0.19     
## ---
## 'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 95%;
## for two-sided hypotheses, the value tested against lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.

Ordinal

df_double$reactance_ordinal <- factor(df_double$reactance,
                                      levels = 0:5, 
                                      ordered = TRUE)

formula <- bf(
  reactance_ordinal ~ 
    persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw +
    
    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +
    pushing_self_cb + pushing_partner_cb +
    
    day +
    
    # Random effects
    (persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw | coupleID),

  autocor = ~ ar(time = day, gr = coupleID:userID, p = 1)
)



prior1 <- c(
  set_prior("normal(0, 2.5)", class = "b"),
  set_prior("normal(0, 5)", class = "Intercept"),
  set_prior("normal(0, 1)", class = "sd", group = "coupleID", lb = 0),
  set_prior("normal(0, 0.075)", class = "ar", lb = -1, ub = 1),
  set_prior("normal(0.5, 2.0)", class = "sderr", lb = 0)
)




#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

reactance_ordinal <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = brms::cumulative(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  iter = 5000,
  warmup = 2000,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", "reactance_ordinal")
)
pp_check(reactance_ordinal, type='hist')
## Using 10 posterior draws for ppc type 'hist' by default.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

pp_check(reactance_ordinal)
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

loo(reactance_ordinal)
## 
## Computed from 12000 by 756 log-likelihood matrix.
## 
##          Estimate   SE
## elpd_loo   -700.9 31.8
## p_loo       162.0 10.5
## looic      1401.9 63.7
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.0, 0.7]).
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. ESS
## (-Inf, 0.7]   (good)     731   96.7%   15      
##    (0.7, 1]   (bad)       23    3.0%   <NA>    
##    (1, Inf)   (very bad)   2    0.3%   <NA>    
## See help('pareto-k-diagnostic') for details.
plot(reactance_ordinal, ask = FALSE)

summarize_brms(
  reactance_ordinal, 
  model_rows_fixed = model_rows_fixed_ordinal,
  model_rows_random = model_rows_random_ordinal,
  model_rownames_fixed = model_rownames_fixed_ordinal,
  model_rownames_random = model_rownames_random_ordinal,
  exponentiate = T) %>%
  print_df(rows_to_pack = rows_to_pack_ordinal)
OR l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercepts
Intercept NA NA NA NA NA NA
Intercept[1] 3.64* 1.78 8.61 1.010 658.58 765.06
Intercept[2] 8.44* 3.84 27.47 1.020 302.81 359.10
Intercept[3] 25.19* 10.03 133.27 1.026 191.87 249.54
Intercept[4] 122.91* 40.86 1297.68 1.032 148.93 239.51
Intercept[5] 6151.51* 960.88 341247.23 1.033 117.18 216.26
Within-Person Effects
Daily persuasion experienced 1.10 0.94 1.29 1.009 2682.51 1046.48
Daily persuasion utilized (partner’s view) 1.01 0.84 1.20 1.008 524.47 454.81
Daily pressure experienced 0.93 0.72 1.20 1.003 2711.17 5306.92
Daily pressure utilized (partner’s view) 1.00 0.69 1.42 1.004 7386.06 4975.11
Daily pushing experienced 1.07 0.90 1.31 1.009 692.83 625.00
Daily pushing utilized (partner’s view) 0.93 0.74 1.16 1.003 3195.70 1990.83
Day 0.68 0.29 1.48 1.005 1929.74 2156.18
Daily weartime NA NA NA NA NA NA
Between-Person Effects
Mean persuasion experienced 0.62 0.21 1.75 1.006 775.50 1753.38
Mean persuasion utilized (partner’s view) 1.61 0.55 4.80 1.011 408.84 598.36
Mean pressure experienced 1.72 0.51 5.91 1.007 489.55 314.92
Mean pressure utilized (partner’s view) 0.81 0.23 2.88 1.009 405.01 653.48
Mean pushing experienced 0.66 0.14 3.01 1.001 3194.62 5159.99
Mean pushing utilized (partner’s view) 0.64 0.12 3.26 1.004 1109.14 1145.95
Mean weartime NA NA NA NA NA NA
Random Effects
sd(Intercept) 1.37 0.92 2.18 1.03 153.00 57.31
sd(Daily persuasion experienced) 0.24 0.01 0.57 1.00 1351.72 2078.86
sd(Daily persuasion utilized (partner’s view)) 0.18 0.01 0.48 1.00 2111.36 4767.30
sd(Daily pressure experienced) 0.93 0.51 1.52 1.01 1090.67 2141.89
sd(Daily pressure utilized (partner’s view)) 0.38 0.01 1.18 1.00 2797.62 5301.03
sd(Daily pushing experienced) 0.22 0.01 0.67 1.01 331.05 134.09
sd(Daily pushing utilized (partner’s view)) 0.16 0.01 0.55 1.00 2639.35 2745.31
Additional Parameters
ar[1] 0.00 -0.14 0.15 1.00 6348.79 3182.63
nu NA NA NA NA NA NA
shape NA NA NA NA NA NA
sderr 0.82 0.03 2.65 1.05 106.36 100.04
sigma NA NA NA NA NA NA
disc 1.00 1.00 1.00 NA NA NA 
hypothesis(reactance_ordinal, "pressure_self_cw > pushing_self_cw")
## Hypothesis Tests for class b:
##                 Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio Post.Prob Star
## 1 (pressure_self_cw... > 0    -0.14      0.18    -0.44     0.15       0.28      0.22     
## ---
## 'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 95%;
## for two-sided hypotheses, the value tested against lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.

Binary

introduce_binary_reactance <- function(data) {
  data$is_reactance <- factor(data$reactance > 0, levels = c(FALSE, TRUE), labels = c(0, 1))
  return(data)
}



df_double <- introduce_binary_reactance(df_double)
if (use_mi) {
  for (i in seq_along(implist)) {
    implist[[i]] <- introduce_binary_reactance(implist[[i]])
  }
}


formula <- bf(
  is_reactance ~ 
    persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw +
    
    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +
    pushing_self_cb + pushing_partner_cb +
    
    day +
    
    # Random effects
    (persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw | coupleID),

  autocor = ~ ar(time = day, gr = coupleID:userID, p = 1)
)



prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b"),
  brms::set_prior("normal(0, 10)", class = "Intercept", lb=0, ub=5), # range of the outcome scale
  
  brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0),

  brms::set_prior("cauchy(0, 5)", class = "ar", lb = -1, ub = 1)
  #brms::set_prior("cauchy(0, 10)", class = "sigma", lb = 0)
)


#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

is_reactance <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = bernoulli(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  iter = 5000,
  warmup = 2000,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", "is_reactance")
)
pp_check(is_reactance, type='hist')
## Using 10 posterior draws for ppc type 'hist' by default.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

pp_check(is_reactance)
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

try(loo(is_reactance))
## 
## Computed from 12000 by 756 log-likelihood matrix.
## 
##          Estimate   SE
## elpd_loo   -415.8 15.2
## p_loo       356.0 13.7
## looic       831.5 30.5
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.7, 1.2]).
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. ESS
## (-Inf, 0.7]   (good)       4    0.5%   504     
##    (0.7, 1]   (bad)      381   50.4%   <NA>    
##    (1, Inf)   (very bad) 371   49.1%   <NA>    
## See help('pareto-k-diagnostic') for details.
plot(is_reactance, ask = FALSE)

summarize_brms(
  is_reactance, 
  model_rows_fixed = model_rows_fixed,
  model_rows_random = model_rows_random,
  model_rownames_fixed = model_rownames_fixed,
  model_rownames_random = model_rownames_random,
  exponentiate = T) %>%
  print_df(rows_to_pack = rows_to_pack)
OR l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 1.03 0.16 6.23 1.000 5639.90 7417.58
Within-Person Effects
Daily persuasion experienced 1.42 0.77 2.83 1.001 4607.59 6137.24
Daily persuasion utilized (partner’s view) 1.04 0.47 2.31 1.001 4472.27 6807.48
Daily pressure experienced 1.03 0.33 3.42 1.001 4972.87 6595.39
Daily pressure utilized (partner’s view) 0.87 0.19 3.83 1.000 5064.79 6729.28
Daily pushing experienced 1.33 0.63 3.01 1.001 4555.22 5161.88
Daily pushing utilized (partner’s view) 0.70 0.24 1.78 1.001 4672.92 5808.05
Day 0.47 0.03 7.93 1.000 6355.54 7827.73
Daily weartime NA NA NA NA NA NA
Between-Person Effects
Mean persuasion experienced 0.52 0.03 10.43 1.001 5678.41 7315.27
Mean persuasion utilized (partner’s view) 3.05 0.16 59.43 1.000 5731.55 6804.92
Mean pressure experienced 1.84 0.07 46.64 1.000 6006.56 7426.18
Mean pressure utilized (partner’s view) 1.08 0.04 30.40 1.001 5983.48 7130.34
Mean pushing experienced 0.55 0.01 31.44 1.000 7338.37 8308.91
Mean pushing utilized (partner’s view) 0.33 0.01 18.09 1.000 7538.67 8622.63
Mean weartime NA NA NA NA NA NA
Random Effects
sd(Intercept) 6.22 4.13 8.67 1.00 3353.77 4628.04
sd(Daily persuasion experienced) 1.28 0.17 2.67 1.00 1555.96 1586.65
sd(Daily persuasion utilized (partner’s view)) 1.18 0.07 2.75 1.00 1318.93 2480.79
sd(Daily pressure experienced) 2.98 1.24 5.15 1.00 2559.40 2743.37
sd(Daily pressure utilized (partner’s view)) 1.20 0.05 3.49 1.00 4330.74 5308.94
sd(Daily pushing experienced) 0.81 0.04 2.11 1.00 2226.82 4965.39
sd(Daily pushing utilized (partner’s view)) 0.74 0.03 2.21 1.00 3602.74 4943.96
Additional Parameters
ar[1] 0.11 -0.08 0.31 1.00 1614.73 2746.80
nu NA NA NA NA NA NA
shape NA NA NA NA NA NA
sderr 7.40 4.37 11.33 1.00 1951.96 3229.35
sigma NA NA NA NA NA NA 
hypothesis(is_reactance, "pressure_self_cw > pushing_self_cw")
## Hypothesis Tests for class b:
##                 Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio Post.Prob Star
## 1 (pressure_self_cw... > 0    -0.26      0.74    -1.47     0.96       0.56      0.36     
## ---
## 'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 95%;
## for two-sided hypotheses, the value tested against lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.

Report All Models

if (report_ordinal) {
  model_rows_random_final <- model_rows_random_ordinal
  model_rows_fixed_final <- model_rows_fixed_ordinal
  model_rownames_fixed_final <- model_rownames_fixed_ordinal
  model_rownames_random_final <- model_rownames_random_ordinal
  rows_to_pack_final <- rows_to_pack_ordinal
} else {
  model_rows_random_final <- model_rows_random
  model_rows_fixed_final <- model_rows_fixed
  model_rownames_fixed_final <- model_rownames_fixed
  model_rownames_random_final <- model_rownames_random
  rows_to_pack_final <- rows_to_pack
}



all_models <- report_side_by_side(
  pa_sub,
  pa_obj_log,
  mood,
  reactance,
  is_reactance,
  
  model_rows_random = model_rows_random_final,
  model_rows_fixed = model_rows_fixed_final,
  model_rownames_random = model_rownames_random_final,
  model_rownames_fixed = model_rownames_fixed_final
) 
## [1] "pa_sub"
## [1] "pa_obj_log"
## [1] "mood"
## [1] "reactance"
## [1] "is_reactance"
# pretty printing
summary_all_models <- all_models %>%
  print_df(rows_to_pack = rows_to_pack_final) %>%
  add_header_above(
    c(" ", "Subjective MVPA" = 2, 
      "Device-Based MVPA" = 2, 
      "Mood" = 2,
      "Reactance Gaussian" = 2, 
      "Reactance Dichotome" = 2)
  )

export_xlsx(
  summary_all_models, 
  rows_to_pack = rows_to_pack_final,
  file.path("Output", "AllModels.xlsx"), 
  merge_option = 'both', 
  simplify_2nd_row = TRUE,
  colwidths = c(38, 7.2, 13.3, 7.2, 13.3,7.2, 13.3,7.2, 13.3,7.2, 13.3),
  line_above_rows = c(1,2),
  line_below_rows = c(-1)
)

summary_all_models
Subjective MVPA
Device-Based MVPA
Mood
Reactance Gaussian
Reactance Dichotome
IRR pa_sub 95% CI pa_sub exp(Est.) pa_obj_log 95% CI pa_obj_log b mood 95% CI mood b reactance 95% CI reactance OR is_reactance 95% CI is_reactance
Intercept 26.12* [18.60, 36.84] 117.56* [106.28, 129.91] 4.92* [ 4.71, 5.13] 0.54* [ 0.27, 0.81] 1.03 [0.16, 6.23]
Within-Person Effects
Daily persuasion experienced 0.99 [ 0.89, 1.10] 1.00 [ 0.97, 1.02] 0.02 [-0.01, 0.05] 0.03 [-0.02, 0.08] 1.42 [0.77, 2.83]
Daily persuasion utilized (partner’s view) 0.93 [ 0.84, 1.03] 1.02* [ 1.00, 1.05] -0.02 [-0.05, 0.01] 0.00 [-0.06, 0.06] 1.04 [0.47, 2.31]
Daily pressure experienced 1.02 [ 0.78, 1.39] 0.98 [ 0.93, 1.04] -0.03 [-0.11, 0.04] -0.02 [-0.12, 0.07] 1.03 [0.33, 3.42]
Daily pressure utilized (partner’s view) 0.64* [ 0.51, 0.82] 1.01 [ 0.96, 1.07] 0.04 [-0.04, 0.11] 0.00 [-0.13, 0.12] 0.87 [0.19, 3.83]
Daily pushing experienced 1.03 [ 0.89, 1.22] 1.01 [ 0.98, 1.04] -0.05* [-0.10, -0.01] 0.03 [-0.03, 0.09] 1.33 [0.63, 3.01]
Daily pushing utilized (partner’s view) 1.12 [ 0.96, 1.31] 0.98 [ 0.95, 1.01] 0.01 [-0.04, 0.05] -0.02 [-0.10, 0.06] 0.70 [0.24, 1.78]
Day 1.07 [ 0.78, 1.47] 1.00 [ 0.92, 1.08] -0.08 [-0.22, 0.07] -0.15 [-0.42, 0.12] 0.47 [0.03, 7.93]
Daily weartime NA NA 1.00 [ 1.00, 1.00] NA NA NA NA NA NA
Between-Person Effects
Mean persuasion experienced 1.22 [ 0.73, 2.08] 0.88 [ 0.74, 1.04] 0.27 [-0.04, 0.59] -0.12 [-0.49, 0.26] 0.52 [0.03, 10.43]
Mean persuasion utilized (partner’s view) 1.25 [ 0.73, 2.17] 1.07 [ 0.89, 1.27] 0.08 [-0.23, 0.40] 0.16 [-0.24, 0.56] 3.05 [0.16, 59.43]
Mean pressure experienced 1.79 [ 0.89, 3.65] 1.12 [ 0.90, 1.40] -0.01 [-0.37, 0.35] 0.18 [-0.25, 0.61] 1.84 [0.07, 46.64]
Mean pressure utilized (partner’s view) 0.52 [ 0.25, 1.12] 0.77* [ 0.62, 0.95] -0.11 [-0.51, 0.27] -0.09 [-0.57, 0.40] 1.08 [0.04, 30.40]
Mean pushing experienced 0.40* [ 0.16, 1.00] 1.09 [ 0.82, 1.45] -0.46 [-0.98, 0.07] -0.17 [-0.76, 0.42] 0.55 [0.01, 31.44]
Mean pushing utilized (partner’s view) 0.57 [ 0.26, 1.29] 1.23 [ 0.97, 1.57] -0.19 [-0.65, 0.26] -0.11 [-0.71, 0.52] 0.33 [0.01, 18.09]
Mean weartime NA NA 1.00* [ 1.00, 1.00] NA NA NA NA NA NA
Random Effects
sd(Intercept) 0.67 [ 0.49, 0.89] 0.27 [0.20, 0.35] 0.63 [0.49, 0.81] 0.58 [ 0.44, 0.77] 6.22 [ 4.13, 8.67]
sd(Daily persuasion experienced) 0.20 [ 0.05, 0.35] 0.06 [0.03, 0.09] 0.03 [0.00, 0.07] 0.06 [ 0.00, 0.15] 1.28 [ 0.17, 2.67]
sd(Daily persuasion utilized (partner’s view)) 0.19 [ 0.04, 0.34] 0.05 [0.02, 0.08] 0.05 [0.00, 0.11] 0.04 [ 0.00, 0.12] 1.18 [ 0.07, 2.75]
sd(Daily pressure experienced) 0.16 [ 0.01, 0.48] 0.06 [0.00, 0.15] 0.09 [0.00, 0.26] 0.43 [ 0.27, 0.66] 2.98 [ 1.24, 5.15]
sd(Daily pressure utilized (partner’s view)) 0.14 [ 0.01, 0.43] 0.04 [0.00, 0.10] 0.15 [0.01, 0.33] 0.21 [ 0.01, 0.58] 1.20 [ 0.05, 3.49]
sd(Daily pushing experienced) 0.24 [ 0.02, 0.52] 0.08 [0.01, 0.15] 0.09 [0.01, 0.16] 0.09 [ 0.00, 0.22] 0.81 [ 0.04, 2.11]
sd(Daily pushing utilized (partner’s view)) 0.15 [ 0.01, 0.35] 0.04 [0.00, 0.10] 0.07 [0.01, 0.16] 0.04 [ 0.00, 0.14] 0.74 [ 0.03, 2.21]
Additional Parameters
ar[1] 0.02 [-0.94, 0.94] 0.28 [0.24, 0.31] 0.45 [0.42, 0.48] 0.01 [-0.07, 0.10] 0.11 [-0.08, 0.31]
nu NA NA NA NA NA NA NA NA NA NA
shape 0.14 [ 0.13, 0.14] NA NA NA NA NA NA NA NA
sderr 0.05 [ 0.00, 0.13] NA NA NA NA NA NA 7.40 [ 4.37, 11.33]
sigma NA NA 0.55 [0.54, 0.57] 0.87 [0.85, 0.89] 0.94 [ 0.89, 1.00] NA NA 
report::report_system()

Analyses were conducted using the R Statistical language (version 4.4.1; R Core Team, 2024) on Windows 11 x64 (build 22635)

report::cite_packages()